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Abstract 

We consider the phenomenon of collapse in the critical Keller-Segel equation (KS) which models 
chemotactic aggregation of micro-organisms underlying many social activities, e.g. fruiting body de- 
velopment and biofilm formation. Also KS describes the collapse of a gas of self-gravitating Brownian 
particles. We find the fluctuation spectrum around the collapsing family of steady states for these equa- 
tions, which is instrumental in derivation of the critical collapse law. To this end we develop a rigorous 
version of the method of matched asymptotics for the spectral analysis of a class of second order differ- 
ential operators containing the linearized Keller-Segel operators (and as we argue linearized operators 
appearing in nonlinear evolution problems). We explain how the results we obtain are used to derive the 
critical collapse law, as well as for proving its stability. 

Key words: Matched asymptotics, critical Keller-Segel equation, collapse and formation of singularities, 
linearized operators. 



1 Introduction 

Phenomena of blowup and collapse in nonlinear evolution equations are hard to simulate numerically and 
the rigorous theory, or at least a careful analysis, is pertinent here. The recent years witnessed a tremendous 
progress in developing of such theories. We can now describe the shape of blowup profile and contraction 
law in Yang-Mills, a-model, nonlinear Schrodinger and heat equations ([DI2J|s2II103I3IE1ISI]I3])Q Yet, 
after 40 years of intensive research and important progress, we still cannot give a rigorous description of 
collapse in the Keller-Segel equations modeling chemotaxis. (See [TTJ Q21 [13j HH [HI [TBI H3 HB1 OS] f° r some 
recent works, 20 , for a nice discussion of the subject, and [2TJ [22J, [23l EH [25] for reviews.) 
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This is not say that the Keller-Segel equations are harder than Yang-Mills, a— model and nonlinear 
Schrodinger equations, they are not, but neither are they less important. They model chemotaxis, which is 
a directed movement of organisms in response to the concentration gradient of an external chemical signal 
([26]. see also (22)- The chemical signals can come from external sources or they can be produced by 
the organisms themselves. The latter situation leads to aggregation of organisms and to the formation of 
patterns and is the case modeled by the Keller - Segel equations. Chemotaxis is believed to underly many 
social activities of micro-organisms, e.g. social motility, fruiting body development, quorum sensing and 
biofilm formation. A classical example is the dynamics and the aggregation of Escherichia coli colony under 
the starvation condition Another example is the Dictyostelium amoeba , where single cell bacterivores, 
when challenged by adverse conditions, form multicellular structures of ~ 10 5 cells [28j[29]. Also endothelial 
cells of humans react to vascular endothelial growth factor to form blood vessels through aggregation [30J . 

We assume that the organism population is large and the individuals are small relative to both the 
domain, ft C K d (d = 1, 2, 3) as well as the typical distance between the organism is much larger than their 
size. One can derive in the mean-field approximation the Keller-Segel system governing the organism density 
p : il x M + — y R + and chemical concentration c : ft x R + — s- M + J2E1 [3T] • As the chemical diffuses much faster 
than organisms, one makes a simplifying assumption of instantaneous interaction (adiabatic assumption) 
which, after rescaling and a minor simplification, leads the Keller-Segel equations to the form 



with p and c satisfying the no-flux Neumann boundary conditions. The equations ([lj appear also in the 
context of stellar collapse (see [3J2 (33J [34j [35]); similar equations - the Smoluchowski or nonlinear Fokker - 
Planck equations - models non- Newtonian complex fluids (see [Ml EZ1 GUI 132] ■ 

Arguably, the most interesting feature of the Keller-Segel equations is that they can develop, in finite 
time, infinite mass at a point in space. As argued below, the 'collapsing' profile and contraction law have 
a universal (close to self-similar) form, independent of particulars of initial configurations and, to a certain 
degree, of the equations themselves, and can be associated with chemotactic aggregation. Though the 
equations are rather crude and unlikely to produce patterns one observes in nature or experiments, the 
collapse phenomenon could be useful in verifying assumptions about biological mechanismso 

We now concentrate on the (energy) critical case of d = 2 and fi = K 2 . It was shown in [331 EH] that 
solutions of ([T]) with the mass 



blow up in finite time. Ref. [13 exhibited blowup solutions with explicit blowup rate and explicit asymptotics, 
which was confirmed in J45, 46 by a different technique relying on results of the present paper. However, 
the problem of describing the dynamics of blowup, i.e. blowup rate and profile for an open set of initial 
conditions is still open. As is shown below, this paper makes a considerable progress toward its solution. 

Of a critical importance here are the following key properties of the equation <jTJ) : 
• It is invariant under the scaling transformations p(x, t) — > j?p (jx, -p-t) and c(x, t) — > c (jx, -p-i) . 

2 There are numerous refinements of the Keller - Segel equations, e.g. taking into account finite size of organisms (;401 |41II42| ) 
preventing the complete collapse, which model the chemotaxis more precisely. We believe techniques we outline and develop 
here can be applied to these models as well. 




(1) 
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• It has the static solution, R(x) := ( 1+ u i s)5 > C(x) := — 21n(l + \x\ 2 ). 

• The total 'mass' is conserved: J p(x,t)dx — const H 

The stationary solution R(x) has the total mass J R(x)dx = 87r, which is exactly the sharp threshold between 
global existence and singularity development in solutions to ([T]), mentioned above. 

The properties above yield that ([l) has in fact the family of static solutions \~ 2 R(x/\), C(x/X), A > 0@ 
and suggest a likely scenario of collapse: sliding along this family in the direction of A — > 0. Indeed, we 
conjecture that, like in Struwe's result |49] for equivariant wave maps from the Minkovskii space-time, M 2+1 , 
to the 2-sphere, S 2 , for any solution, p(x,t), of JT]), collapsing up at time T, there are sequences Ai — > and 
U — > T, s.t. p(Xiy,ti) converges to the stationary solution R{y), as i — > oo. Thus the most interesting and 
natural initial conditions for (JlJ are those close to the manifold {R\(x)\X > 0}. 

This discussion brings us to the first step of the theory of collapse in the Keller - Segel system - deter- 
mining the low- lying spectrum of fluctuations around the family R\ (x) . This would determine whether this 
family is stable. In this paper we find this spectrum and to do this we develop a rigorous version of the 
method of matched asymptotics. 

Now, we discuss, following [3S], a natural approach to this problem. Since the blowup profile is expected 
to be radially symmetric, it is natural to start with radially symmetric solutions. In this case, the system 
(JT|), which consists of coupled parabolic and elliptic PDEs, is equivalent to a single PDE. Indeed, the change 
of the unknown, by passing from the density, p(x,t), to the normalized mass, 

m{r,t) ■= ^~ I p(x,t) dx, 

l ^ J\x\<r 

of organisms contained in a ball of radius r, discovered by [52j [20], maps two equations (QJ into a single 
equation 

d t m = 3®'m + r~ 1 md r m, (2) 

on (0, oo) (with initial condition m Q (r) := ^ J\ x \< r Po( x ) dx). Here is the n-dimensional radial Laplacian, 
A^™' := r^' 1-1 ) d r r n ~ 1 d r = d 2 + ^r-d r . Thus ([1]) in the radially symmetric case is equivalent to ([2]). 

The equation @ has the following key properties, inherited from the corresponding properties of ([l])): 

• It is invariant under the scaling transformations m(r,i) — > m (jr, j^t) . 

• It has the static solution x( r ) nip ■ 

• The total 'mass' is conserved: 2it\im r ^ 00 m(r,t) — J p{x,t)dx — const. 

3 Another important property of JTJ that it is a gradient system with the (free) energy F(p) = J^i^P A - 1 p + plnp)dx, 
which plays the key role in other papers, is not used in our approach. 

4 It seems this family was discovered in 1471 . It is shown in |48| that belongs to the two parameter family Hi (x) := [r/X), 

where R^(x) := 2(fi — 2) 2 rrjr^ m-ihii > P > 2. Our case is )i = 4. If 2 < fj, < 4, then the mass at the origin is non- 

(1+ 1 x I ) 

zero, and if [i > 4, then the mass at the origin is negative and hence the static solution is not physical. For [i = 4, due 
to the sharp logarithmic Hardy-Littlewood-Sobolev inequality, these static solutions unique and minimize the free energy, 
F(p) = J r2 {\p A _1 p + p lnp) dx, for the fixed mass J p = p ( [501 1511 ). We conjecture that the same is true for 2 < p < 4. 
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As in the case of ([T]), the properties above yield the manifold of static solutions Mo := {x(r/X) | A > 0} 
and suggest a likely scenario of collapse: sliding along Mo in the direction of A — > 0. To analyze the collapse, 
we pass to the reference frame collapsing with the solution, by introducing the adaptive blowup variables, 

r f* 1 

m(r, t) — u(y, r), where y — — and r = / „. . ds, 

X Jo * K s ) 

where A : [0,T) — > [0, oo), T > 0, is a positive differentiable function (compression or dilatation parameter), 
s.t. X(t) — > as t j" T . The advantage of moving to blowup variables is that the function u is expected 
to have bounded derivatives and the blowup time is eliminated from consideration (it is mapped to oo). 
Writing ([2|) in blowup variables, we find the equation for the rescaled mass function 

d T u = A^'u + y~ 1 ud y ii — ayd y u, (3) 

where a := —XX. 

To investigate stability properties of the rescaled stationary solution \(y) , we decompose solutions u(y, r) 
of equation ([3]) into the main term, x(y)-> an d the fluctuation </>(y,r), u(y,r) — x{y) + 4>{v,t). Substituting 
this decomposition into ([3]) gives the equation for the fluctuation 0, 

d T cj) = -L a cb + F a + N(cb), (4) 

where the forcing and nonlinear terms are F a :— —j^i^yi and N((fi) :— ^(f>d y (f>, and the linear operator, 
L a is given by 

^ ■= -A^ - ir A_ + ^-^dy + ayd y . (5) 

An important fact here is that the operator L a is self-adjoint on the space L 2 (M. + ,'y a (y)y 3 dy), where 
la 1 ^ 2 (y) = x{y) e ^ v i with the inner product (f,g) :~ f °° f(y)g(y) Ja(y)y 3 dy. One can check the self- 
adjointness of L a directly or use the unitary map 

<Kv) -> 7a /2 (y)^(2/), (6) 

from L 2 ([0, oo), r y a (y)y 3 dy) to L 2 ([0, oo), y 3 dy), which maps this operator L a into the operator £ a := 
lV 2 L a la 1 ^ 2 , acting on the space £ 2 ([0, oo), y 3 dy). The latter operator can be explicitly computed to be 

L a := -AW - — « + W + j^-2a. (7) 
(1 + y z y 4 L + y 

This operator is of the Schrddinger type with the real continuous potential tending to oo as y —> oo. 
Therefore, by standard arguments (see e.g. [S3]), it is self-adjoint and its spectrum is purely discrete. Hence 
L a is self-adjoint on the space L 2 ([0, oo), j a (y)y 3 dy) and has purely discrete spectrum as well. Going through 
with our analysis shows that a(r) — > as r — > oo, which actually complicates the problem and which tells us 
that the collapse is slower than parabolic one, X(t) = ^/ao(T — t), for which a(r) = —X(t)X(t) is a constant 
(say, ao). 

Now, it is clear that the stability of the the profile x(v) is determined largely by the spectrum of the 
operator C a . If the operator C a has strictly positive spectrum, one expect the fluctuations <f> will die out 
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as t — y oo and consequently the solution of ([3]) will tend to x(j/)i while the solution of will approach 
x(r/X{t j). On the other hand, if the operator C a has negative eigenvalues then one expects instability. The 
latter though is always the case, since the equations have a negative scaling mode (for a fixed parabolic 
scaling it is connected to possible variation of the blowup time.) □ 

If the number of negative eigenvalues is finite, say k, then one either goes to an invariant manifold theory 
and constructs the central-unstable manifold or, uses the (related) modulation theory and embeds x(r/A) 
into a k— parameter family of almost solutions, say Xp( r M)? where p stands for the k — 1 parameters (with 
A, or a, counted as the first parameter), chosen so that the tangent space of the deformation (or almost 
center-unstable) manifold M. :— {Xp( r M) I A > 0, p} at Xpiu) is equal approximately to the eigenspace of 
negative and (almost) zero spectrum of C a . Then we can choose the parameters p = p(r) and a — a(r) (or 
A = A(t)), so that the solution u(y,r) can be decomposed as 

u{y,r) = X P (r)(y) + <P(v,t), (8) 

with the fluctuation </>(y,r) orthogonal to the tangent space of M. at X P (y), (dpXp(r)(:), <!>(■, T )) — 0, and 
therefore (approximately) orthogonal to the negative and almost zero spectrum eigenfunctions of C a - If 
we find such a deformation, Xp(t){u)i then the stability is restored and the solution to (O approaches this 
family as r — > oo. The latter is a big if and this is where the understanding the negative and almost zero 
spectrum eigenfunctions of C a helps. (One should keep in in mind that the deformation of xiu) wu l change 
the linearized operator C a and the gauge transformation ([6]), but both can be easily handled.) 

Theorem[T]below implies that the operator C a of the Keller-Segel system has one negative (corresponding 
to the scaling mode mentioned above) and one near zero eigenvalue, while the third eigenvalue, 2a + + 

O (a In -2 — ) , is positive, but vanishing as a — ¥ 0. (It also isolates the correct perturbation (adiabatic) 
parameter - j^-r-) Hence we have to construct a one-parameter deformation of x(y) (remember that A, or 

a, is counted as a parameter). For technical reasons it is convenient to use a two-parameter family, Xbc{y), 
with an extra relation between the parameters a, b and c. In [48] we choose the family 

Xbc(y) ■= -^7-3. (9) 

with b > 1 and both parameters b and c are close to 1. Note that this family evolves on a different spatial 
scale than cf>(y,T) in ([5]). as it can rewritten as Xbc(y) = Xi. The contraction law is obtained by using 

the orthogonality condition, (dhcXba 4>) = 0. The latter is equivalent to two conditions, 

d T (dbcXb(r)c(r){-), </>(', T)) =0 (10) 

and (<9fc c Xb(-r)c(-r)(')j 4>{'i T )) |t=o = 0. We evaluate (jXOj) by using the evolution equation, d T 4> = —L a b c <f> + 
F a bc + N ((/)), for </>, similar to (j4]), which follows by plugging the decomposition (jHJ into ([3]), and the explicit 
expressions, 



1 y 2 1 y 2 

Cbd(y) ■= -rdbXbciy) = , n , Cbc2(y) ■= -rrdcXbdy) = -7— — (11) 
4 c + y z 4o (c + y ) 



5 A similar analysis applies also in the subcritical case M < 8n where the solution converges to a self-similar one as r — > 00, 
which vanishes as t — > 00. In this case the operator C a has strictly positive spectrum. 
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for the tangent vectors to the manifold M.. (The scaling mode 0>co(j/) '■= gj^ydyXbc(y) = r^jjrp a multiple of 
Cbcsiy) which confirms that one of the parameters is superfluous.) This gives ordinary differential equations 
for a, b and c with higher order terms depending on <f>: 



c T + S((j>, a, b, c)a T = 2a — j^xy + ^(07 a : 6, c), 
£ - a, b, c)a T = - rMr + R(cj>, a, b, c), 



(12) 

where d:= 6-1, \S(<t>,a,b,c)\ < J + ?^ } and \R(4>,a,b,c)\ < I -^ yi _^ y [d||0|| L2 + ||(1 + t/ 2 )- These 

higher order terms are controlled by using a differential inequality for the Lyapunov functional 4> i— >■ ||0|| 2 2 
and the inequality (</>,£ a b c 0) > 2a||(/>|| 2 , which follows from our result below. @ (One can also use higher 
order Lyapunov functionals like ((j>,C^ bc (f>), k > 1. The fact that the positive eigenvalues of L a 6 C vanish 
makes estimating (f> a delicate matter.) Finally, we choose a relation between a, b and c so as to eliminate 
large terms in the corresponding vector fields, namely, d — ialn(i), This leads, in the leading order, to the 
differential equation 

(13) 



Hi 

whose solutions, in the leading order, are i (In i + 0(1)) = 2r which results in In = ln2r — lnln2r + 

T^T + ° (eW)- Recalling that \{t)\{t) = -o(r) and using that A(i)A(f) = A(r(i)) _1 9 r A(r(i)), we obtain 
that — In A(r) = J r a(r)dr giving 

, WN (ln2r) 2 (ln2T)lnln2T ^„ 
-lnA(t) = ^- i ^--^ ^ +0(k2r), (14) 

while r is related to £ by 

A(r') 2 dr' = (1/2) [(ln2r)" 1 + 0((ln2r)- 2 lnln2T)] 
x exp 



- ^ ln2T ' )2 + (In2r)(lnln2r) + 0(ln2r) 



(15) 



where the integral was computed asymptotically in a limit t» 1 and T is a constant of integration determined 
by initial conditions. Solving the equations (fT4")l and (fl5|) together for t — > T yields the law 



A(i) = (T-^e-l^-'^fcx+oCl)), (16) 

which coincides in the leading order (up to the constant Ci) with the one obtained in [451 125] ■ The 
constant c\ can be obtained only if we consider a next order correction beyond the accuracy of the equation 
(fT5| (see HB]) which is outside the scope of this paper. 

The above arguments show that the spectral analysis of the linearized equation on the the collapse or 
blowup profile is the key step in describing critical collapse or blowup laws for nonlinear evolution equations. 
(This also applies to stability analysis of stationary and traveling wave solutions.) Typically, this is a rather 



^Presently, without taking into account the nonlinearity in the equation d T <f> = —L a i, c <f> + F a b c + N(<j>) ( 48 ). 
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subtle affair with very few general techniques available. In this paper we develop such a technique for 
differential operators of the form (J7J , 



C 



-AW 



1 



(i + y 2 ) 2 



2 2 

a y 



W a {y), 



(17) 



defined on the space L 2 ([Q, oo), y 3 dy), which is the subspace of radially symmetric functions in L 2 (R 4 ). Here, 
as indicated above, the parameter a is assumed small and positive and the potential W a , to satisfy the bound 



0< W a (x) < 



Ca 
1 + y 2 



(18) 



where C is a positive constant. We assume that W a (x) is positive in order to fix the bottom of the spectrum: 
now C > (see below). Our main result is the derivation of an approximate equation for the low- lying 
eigenvalues of C, which enter into the stability analysis mentioned above. Let ^Sf be the digamma function, 
defined by Vt'(s) = 4- lnT(s), where T(s) is the gamma function. Define 



and K := In 2 



< fi 

-m) 



W a (y) 

,2\2 



y dy < Ca. 



/o (1 + ?/) 2 

0.577216 ... is the Euler-Mascheroni constant. We have 



1 — 27, where 7 

Theorem 1. The operator C is self-adjoint on L 2 ([0, 00), y 3 dy), positive and its spectrum is discrete. For a 
small, the eigenvalues of C in the interval [0,Ca], for any given C > 0, satisfy the equation 



A 



mi 

a 



* 1 



1 + ( al/2ln ^) 



(19) 



As a — > 0, solutions of 



A 



In 



1 



* 1 



X_ 

2a 



K 



1 



converge to the eigenvalues of C We solve this equation approximately in Section [5] to obtain 

n = 



A n — 



, ttt, +0(aln~ 3 i 
2na 



2d 



In ■ 



-_R-+ 7 -ff„_i- 



7HTT 



Ofaln" 3 i) n>l, 



(20) 



(21) 



where H n := Ylk=i V s - These approximations to the eigenvalues, especially the one obtained by solving 
numerically (|20p. match remarkably well with the numerical computation of the spectrum of C, the results 
of which are given in Figure [T] below. The fact that the numerical solution to (|2Tfl) gives much better 
approximation to the true eigenvalues is not surprising: the approximation (1211) has the logarithmic error 
while the equation (|20|) is obtained with the power accuracy. 



Figure Q] compares the eigenvalue approximations obtained using (|20[) and (|21l) against the numerical 
computation of the first three eigenvalues of C. We have taken W a = 2a/(l + y 2 ) (this gives fi = a). 
Numerical procedure is explained in Section [6] The computations confirm the spectral picture we have 
obtained analytically with the high precision numerical computations. 

We analyze the spectrum of C in the interval [0, Ca] for any fixed C independent of a. This is sufficient 
for the stability analysis for the problem described above. However, we believe that our results are valid in 
a larger interval. 



KSspec, September 21, 2011 



8 



A/a 

5 



1. xlO" 



1. xlO" 



O.OOOOl 



0.0001 



0.001 



0.01 



Figure 1: The quantity Aj/a is plotted against a for the first three eigenvalues obtained in three different 
ways. The solid lines are numerically computed, the dashed lines are obtained using the expressions in 
(|21[) . and the dotted line (coinciding with the numerical computation for the first eigenvalue) is obtained by 
numerically solving the eigenvalue equation (l/TUl) . The eigenvalues plots obtained by solving the full equation 
and the equation (|20p are very close. 



We have already discussed the self-adjointness of C on L 2 ([0, oo), y 3 dy) and the discreteness of its spec- 
trum. The scaling symmetry of ([1]) implies that the function rji(y) := 1/(1 + y 2 ) is a null vector of the 
operator 

Co := -A^ - " 



1 



y 



2\2 ' 



By the Perron-Frobenius argument this implies that the above operator is non-negative and has the non- 
degenerate eigenvalue at 0. The first fact implies that C > 0. 

As was mentioned above, TheoremQ]is proven by a rigorous version of the method of matched asymptotics 
(see [54]). Though this method is standard; other instances of its use in spectral problems can be found in 
[55l l56l l57l l58l [59] . we, however, believe that our extension of this method is novel and robust and can 
be used a large variety of linear differential operators arising in the linearization of nonlinear equations and 
hopefully can be extended to nonlinear ones as well (in this case perturbation series below should be replaced 
by fixed point iterations). 

Indeed, like arguments outlined above, our approach is fairly general. It requires essentially only the 
properties listed above: the scaling symmetry and existence of a stationary solution. (In case of asymptotic 
motion of solitons, the scaling symmetry is replaced by translational, or more generally Galilean or Poicare, 
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invariance.) Indeed, the potential term — rp^sp comes from linearlizing the nonlinear part of the equation 
on the stationary solution x{v) ( m^rp = ^yX(j/))? while the confining potential ^a 2 ?/ 2 comes from the 
term (vector- field), ayd y , generated by the time-dependent rescaling, and it occurs in all such problems. 
We do not use the explicit form of — ( 1+ ^2)2 (besides of its asymptotics at y — > oo and at y — > 0), but 
the fact that, since the stationary solution, x{v)i breaks the scaling symmetry, it leads to the zero mode 
r)i{y) ■— ydyx(y) of the original linearization, Co := 7q //2 £o7o ^ 2 = — ~ [r+^^p (after the transformation 

Finally, we mention the major limitation of our set-up - we deal with radially symmetric solutions. Since 
the only stationary solution is radially symmetric (unfortunately, in contrast to biological observations) , the 
linearized operator for the full equation commutes with rotations and therefore can be decomposed in the 
direct sum of radial operators. Hopefully, our analysis can be extended to each component of this sum. 

This paper is organized as follows. In Section [51 using perturbation theory, we solve the eigenvalue 
problem 

Cc^x = \<j>\ (22) 

in the inner region and then proceed to find the leading order expression. We also use perturbation theory 
in Section [3] to solve (f22j) in the outer region and find the leading order behaviour of this solution. In 
Section [4] we match the inner and outer solutions and in Section [5] we solve the equation ([20)) to obtain the 
solutions (|21[) . Finally, in Section [6] we briefly discuss our numerical procedure. In Appendix A we give 
explicit derivations of some of the expressions of Section [31 which were obtained with reference to the theory 
of special functions. 

In what follows, we use the notation / < g for /, g > 0, if there exists a positive constant C such that 
/ < Cg holds. If the inequality |/| < C\g\ holds then we write / = O (g). We also write / <C g or / = o (g) 
if f / g — > as a or y approach some limit (always specified) and / ~ g if the quotient converges to 1. 



2 Solutions in the Inner Region 

In what follows, A > is a spectral parameter and < a « 1. To simplify the expression we assume in what 
follows that A < Ca for some C > 0. Below we take Ri = £ijah with a <C e, <C ■\fj^fj 
this section is the following 

Proposition 2. If a is small enough, then the solution to the eigenvalue problem 
[0, Ri] is unique, modulo an overall constant factor. For y € [R ,Rj\, R Q 3> 1 as a 
given by 

' ^Alny 2 + \(2\ + n) +K l , 



The main result of 



in the inner region 
► 0, this solution is 



whe 



iii = o a 



y 



In 2 2/ , 1 



— + a 2 y 2 + ef + 
V 



A 2 . 1 
In - 

a a 



-1 2 



(23) 
(24) 



and the O (/) signifies the bound |0 (/) | < C f with a uniform constant C . 
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Proof. For Ri > 0, let Ij-j^ be the norm defined on measurable functions / : [0, E4] — > R by 

11/11,:= sup \(l + y 2 )f\ 

v<R^ 

and let Xi := {/ : [0, Ri] — > K | H/l^ < 00} be the associated Banach space. 

Equation (|22|) can be written as (Ci + Vi)4>\ = 0, where Ci and Vi are defined as 

d 2 3 d 8 Jlr . 1 2 2 „, , , 

4 - - - (TT^I -d Vi := -A + -aV + W a (y). 

The operator is self-adjoint on L 2 ([0, 00), y 3 dy) with essential spectrum cr e;js (£ i ) = [0,oo). It is straight- 
forward to check that the function 

*(») : =TT7 

is a solution to the zero mode equation £,77 = (as mentioned in the introduction, this is due to the scaling 
symmetry of ([!])). Moreover, positivity of 771 and the Perron- Frobenius theorem imply that = inf a(Ci), 
and hence Ci is a positive operator, but is not an eigenvalue. The function 771 is not an eigenvector of Ci 
since it does not lie in L 2 ([0, 00), y 3 dy); we call it a resonance eigenvector. Using that the Wronskian of Ci 
is 1/y 3 , we find a second solution of dr/ = 0: 

m ■= m Qj/ 2 + lny 2 - i?T 2 

This vector is independent of 771, also does not lie in L 2 ([0, 00), y 3 dy) and has a singularity at y = 0. Note 
that i] £ Xi. 

Using variation of parameters we find that the general solution to £j£ = / is = (C~ l f)(y) + 
CiWiiv) + C2^2(y), where Ci and C2 are arbitrary constants and 

(A" 1 /)(y):='7i(y) [ V2(x)f(x)x 3 dx-r l2 (y) f m {x)f(x) x 3 dx. (25) 
Jo Jo 

Lemma 3. Say i/iai a <C £; <C 1- For a small enough, Equation (|22p «as a unique, modulo an overall 
constant factor, solution on [0, Ri] of the form 4>" 1 = 771 + £ wrf/i £ € and 

00 

c = £(-/:r 1 K) n »ft. (26) 

n=l 

T7ie convergence is in Xi and \\C~ 1 Vi\\x^x i ^5 (4 + l) e f m q- 

Remark 1. In fact, we can show convergence of the series in an appropriate norm without the condition on 
the parameter a and the range of y. 

Proof. Recall that C — Ci + Vi. Substituting 0™ = 771 + £ into (|22[) and using £^771 = and the general 
solution of Cit; — f found above, we obtain that £ = — £^~ + Vi-q) + C\i^\ + 6*2772 • The term Ci?7i 
leads to change of an overall factor multiplying 0™ and therefore it can be dropped. Next, if £ G Xi, then 
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so is C i 1 Vi£, (see below). Since 772 ^ X i7 we must, therefore, take C2 = 0, otherwise we would have a 
contradiction. Finally, we rearrange the resulting equation to obtain that 



(I + CrtyiK = -Cr l V iVl . 



(27) 



As (|25|) shows, (|27[) is a Volterra equation and therefore the operator on the r.h.s. has an inverse and this 
inverse can be expanded in the standard perturbation series. We can invert the operator on the left hand 
side on the space X i: provided \\C^ Vi\\xi-yXt < 1- To show the latter property, we compute that 



\C^ Vi fl < ( sup \P{V)^V)\[^§^Z^ Z 



v<Ri 

sup \p(y)m(y)\ 

v<Ri 



y \m{z)v&)\ 3 
P(z) 



(28) 



where p(y) = (1 + y 2 ) and / S Xi. Substituting the expressions for p, rjx, 772 and Vi into the first term and 
using flT5|) gives that 



\PVi\ 



1 L (z) z dz - 



1 



(1 + z 2 ) 2 



1 2 , 2 1 1 
-z 2 + lnz 2 ^ 

2 2z 2 



-A+ T z 2 + Ty Q (z) 



z 3 <iz 



< 



A + a 2 z 2 + 



1 + z 2 



: dz. 



This gives that |p??i| 1 ^ ii (2) z 3 dz < a 2 y 4 + aln(l + y 2 ) + Ay 2 , and hence 

sup \prn\ f '^Q(z)z 3 dz <a 2 R A i +a\nR 2 + XR 2 . 

[0,Ri] Jo P 



(29) 



(30) 



Similarly, we compute that 



and hence, 



iw21 r z " dz ~ ( y2 + (a/ + V(i + in(i + y2)) + a2y6) (i + y2) 

sup |p?72| f V -^-^(z) z 3 dz < a 2 i?f + XR 2 lnifc + (a + A)i? 2 . 

[0,Ri] JO P 



(31) 



Substituting the definition Ri := Ei/a 1 / 2 into (|30[) and (|3Tj) . using 1 3> £j 3> a to simplify In ^- to In - and 
then using the results in (|28|) gives that || ZI" 1 "V^ || js s : i _ i .js s : i < £ 2 ^ In A. and hence a can be taken small enough 
so that HiZ" 1 ViHxi-s-X; < 1- Now inverting the operator on the l.h.s. of (|27p and expanding the inverse into 
the Neumann series completes the proof. □ 

The expression (l23l) for </>\™ is obtained as follows. Due to Lemma|3]and since ||?7i|L = O (1), if a is small 
enough (and £j <C 1), then 



6 A = m - c^Vitn + Ox, 4 + 



-ef In - 

a a 



(32) 
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where the Ox { (e) stands for a function bounded as \\Ox f ( e ) ||x< ^ e - Since |/| < ^r^H/Hi, the latter 



is the same as \Ox { (e) I < e(l + y 2 )" 1 - Substituting the large y expansions rji 



O (y 4 ) and 



rj2 = 1/2 + O (ln(y)/y 2 ) into the expression for C i ViTji (see (|25"]0 and keeping only leading terms needed 
in the region [i? D , gives the expression 



2X + fj, A 



ln ? / + A 



ln 2 y 



iny 



y 2 y 2 



2 2 

a y 



4 4 

Substituting this expression into Q32p. with 771 replaced with its large y expansion and using that 



10 



" A 2, 1" 




'A 2l 1" 

-£i In - 


1 


a a 




a a 





gives ([23]) 



□ 



3 Solutions in the Outer Region 



In the following discussion we take R := e /a^ with a i <C e <C 1 as a — > 0. The main result of this 

in* a 



o3 

ion we tci-tvc -ii .— t /u-^ wiLii - 

1 

section is the following 

Proposition 4. On [i? ,oo) (|22|) /ias a unique solution 0™ , which, for y 6 [i? ,i?j], tofces the form 



1 out 



1 . _ _ _ A A\41 a 



7^ = O I ay 2 | ln(ay 2 )| + ^^(ay 2 ) A"M . (3-1 1 



= In 2/ -In - -In 2 - 1 + 2 7 + * 1 -—--— + ~ + ft , (33) 
a \ 2a J X y z A 

where 

i? 2 

Proof. For a, i? Q and on measurable functions / : [i? D , 00) — > R, we define the norm 

||/|| o := sup ay 2 (ay 2 + l)-^ety 2 / 

Let X a be the corresponding Banach space of functions defined on [R Ol 0) with finite norm ||/|| . 
We write (J22J as (£ + V )<j)\ = 0, where 

C := --^2 - ~ + K(ay) - A and V := U(y) + W a (y). (35) 

«y 2 y dy 

(Recall that V{ay) = \o?y 2 and U(y) = — ■) In the outer region y > we treat V as a small 

potential. The operator C is self-adjoint on L 2 ([0, 00), y 3 dy). 

In what follows the relation / ~ g as y — > 00 means that //g converges to a constant (which might 
depend only on ^-) as y — > 00 and the notation / w g as y — > means that //g converges to 1 as y — > 0. 
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We show in Appendix, Proposition [71 that the equation C <j) = which is the eigenvalue equation for 
the spherically symmetric harmonic oscillator in D = 4, has two linearly independent solutions, <f>o and </>i, 
satisfying the estimates 

A * (ay 2 )' 1 (ay 2 + 1)^,6-^' (36) 



I0o I 



< 



and 



and having the Wronskian 



< 



r i 



A.) 

2d ) 



-(ay 2 + iy 



-i e fv' 



W(<h,<l>i) 



Ay 3 



(37) 



(38) 



Hence 0o £ X a . Using variation of parameters and the Wronskian (|38j) . we find the general solution to 
£ £ = / as £ = C~ l f + C\(f)Q + C-2<f>\, where C± and C2 are arbitrary constants and 



A /*oo />oo 

C- 1 f:=--(^ J fafifdy-faj O 



Lemma 5. /Say </iai ln(£g)| « 1. If the parameter a is small enough, the equation (1221) has a unique, 
modulo and overall constant factor, solution on [R a , 00) of the form 4>q + £, with £ £ X Q and 



(39) 



TTie series converges absolutely in X Q and 



\C- l V \\ Xo ^ Xo <^\\n(aR 2 o )\. 



(40) 



Proof. Substituting </>^ nt = O + £ into (£ + V )<j> = and using that C a (f>o = 0, we obtain £ Q £ + V Q £ = 
— V (j) . Now, using the form of the general solution to £ D £ = / found above and that <\>\ $ X Q gives that 
(I + £~ 1 V )l; = —£~ 1 V (j)o- We note that the choice of the constants C\ = and C 2 = follows from similar 
arguments as in Proposition[3]). The operator on the left hand side can be inverted using the Neumann series 
if ||£~ 1 T4||x ->jc ( , < 1, and hence we estimate LC" 1 ^/ for / e X a : 



K'VJll < - sup | P O | \ 

» y>R a L Jy 



4>lV 



Z 3 dz + 1/30! I 



<foV 



z 3 dz\\\f\\ . 



(41) 



where p(y) :— ay 2 (ay 2 + 1) £a.ei y . For y > R and R large, we have, using (fT8|) and a2_R D <c 1, that 



Using this estimate and (|36l) and (1571) in (14"TT) yields 



zdz + p Vie 1 



Po e S * «k f> ll/llo- 



(42) 
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where po — Po(^f-) and pi = pi(^-) are the prefactors on the r.h.s. of (|36| and (|37|) . respectively. Changing 
the variables of integration as t = and using that aR 2 = e 2 <C 1, we obtain (|40|) . 

Hence, for a is small enough so that the r.h.s. of (14T)|) < 1, we have that ||£~ 1 V^||x ->x < 1, and 
therefore the series converges absolutely, which completes the proof. □ 

Lemma [5] shows that in [i? ,oo), 4>o solves (|22|) in leading order in a: 

rx ut = 4>o + o x (-^\HaRl)\Y (43) 

Next, we show in Appendix that 4>o has the following asymptotics 

A, = ln(«y 2 ) -ln2-l + 2 7 + *fl-^J-^ + | + (ay 2 \n(ay 2 )) , (44) 

where, recall, *S? is the digamma function and 7 = — ^(1) = 0.577216 ... is the Euler-Mascheroni constant. 
Expressions O and (gU) and the observation that I /| < (aj/ 2 ) _1 (ay 2 + l)^e~f y2 \\f\\ Q give (03]) -(EE}- □ 

4 Matching and Eigenvalues 

In this section we prove the main result stated in the introduction. We have two expressions, the inner and 
outer solutions, (|23j) and (|33j) . which solve C<j)\ = \<j)\ in the common region [Ro,Ri\. The inner and outer 
solutions, (1231) and (|33l) . should be equal, up to a constant multiple, in the common region [Rq, R4]. Hence 
we require that 

- ^Alny 2 + \ + -(2\ + ^)+'R l 
4 y 1 4 

= C ( lny 2 - ~ - In - - In 2 - 1 + 2 7 + ~ + * f 1 - + TZ a ) (45) 
\ \y z a A \ 2a J J 

for y e [i? ,i?i]. Here and TZ given in {HJ and (jMj) . 

Note that the remainders in (j2"3")l and (j3"3")l are much smaller than the corresponding leading terms, i. e. 
< a and |ft | < 1, if 

1 



y > In -J]xiR eiRi and ay 2 < 1, (46) 
a 

respectively. We assume that a — > and that ([46]) holds. Equating the leading terms in Eqn (|45|) . i.e. the 
terms which multiples of In y 2 , gives the equation 

C=-^+R c , 

with \Rc\ < m f , where inf is taken over R a < y < Ri satisfying the condition (|46)) . and therefore 

R c = O (a 3 / 2 ). (Here and below we take e% ~ e Q ~ a 1 / 4 .) Similarly, equating the constant terms in (|45|) 
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and substituting the above expression for C gives that 



2A + n = A [ In - +ln2 + 1 - 2 7 - - - * ( 1 - — 
a A V 2a 



R, 



where the higher order term R is bounded as |i?| < inf(a|7£ | + and therefore satishes 

R = o(a^ 2 \n(- 



Rearranging the above equation, assuming that /i + a^O, gives that 



A 



M + a 



, 1 f ^ 

In * 1 

a V 2a 



K 



R 



H + a 



where, recall, K := In 2 — 1 — 27. This is the equation ([T9")) . This proves Theorem [T] 



□ 



5 Solution of fl2DD 

Proposition 6. TTie sei 0/ solutions to (|20l) as a — > is {A n }^L , where X n is given by (I21[) . 

Proof. The term In i on the left hand side of ([20]) is unbounded as a — > whereas the right hand side is 
bounded. Thus, there are two possibilities: either X/(fi + a) <C 1 as a^Oor |A/(a + u)| > C > and there 
is cancelation between In - and ^(1 — A/ 2a). 

We begin with the first case. If X/(fi + a) <C 1 as a — > 0, then, since fi < a, X/2a ~ A/(/x + a) as a -> 0. 
We use this and the fact that *(1) = -7 and *(1 + 5) = *(1) + O (8) to write ((20]) as 



it + a 



In- 



A" + 7 + O 



/ VaA 



(47) 



This equation immediately gives the rough estimate that A/(/x + a) = O (in 1 i) . Substituting this estimate 
into the O (•) term in (|4"?| . then solving the resulting equation for A gives that 



A = 



jU + a 



lni+A + 71 + O (in" 2 i) ' 
Further simplification of the right hand side gives the n = expression in (|21[) . 

If I A/ (// + a) I > C > 0, then there must be cancelation between In ^ and \& (l — in (f2"U|) as already 
mentioned above. The digamma function ty(x) has poles at x = —n for integers n > and hence, for 
cancelation to occur, A/2a must have the form 



A_ 

2a 



l + n + 5. 



where 5 -C 1 as a — > 0. Substituting this form of A/2a into (fT9")l gives the equation 

(1 + n + S) (hx - - m-n - 6) + k] = ^±^. 

V a J 2a 



(48) 



(49) 
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We extract the singular behaviour of n — S) using the identity 

n - 

a , ( _ n _ e) = v [ , ( l_ ( j) + ^ . (50) 

fc=0 

If fc > 1 and <5 < 1/2, then l/(fc + <5) < 1/k + S/k 2 . Using this bound and the fact that = -7 + O (5) 

in J50]) yields that 

1 ™ 1 

*(-n-<5) - 7 + £t + 0(<5). 

fe=i 

Substituting the right hand side for n + S) in (|49| gives the equation 

H + a 



(1 + n + S) I \n--\ + K + j-Y 7 + 0(5) ] = 

V a 6 ti k J 



2a 

As before, a rough estimate of S is In" 1 -. Using this to invert (1 + n + S) gives that 
In 



11^ 1 „ /. _i M U + a 



;=1 fc V a/ 2a(l + n) 1 + O (in" 1 ±) 
and hence, solving this equation for 5, we find that 

S = -L i + o( In' 3 - 



in a + i\ + 7 2^fc=l fc 2a(n+l) 

Substituting this expression into (j48j) gives the eigenvalue approximation f|2 1 [) for n > 1 (with n replaced by 
rt — 1 above) and completes the proof. □ 



6 Numerical Calculation of Spectrum 

To determine eigenvalues and eigenfunctions numerically we used a version of shooting method: we numer- 
ically solved the eigenvalue equation (I2TJ1) with initial conditions (j>\(y = 0) = 1 and -^4>x(y = 0) = for 
each value of A. For general A, solution at ay 2 1 is a linear combination (|54j) . which grows exponentially 
as given by (1571) . We used Newton's method to find values A for which c\ — (i.e. removing exponentially 
growing terms at infinity) in ([54|) . Stopping criterion for Newton's method was to have both <t>\(y) and 
-JU'Pxiy) to decay exponentially for large y. Note that vanishing of 4>\(y) for y -4 00 is not sufficient because 
it would not exclude spurious solution when ci<fio(y) +C2<f>i(y) = at one point only. We controlled numerical 
precision by the matching of numerical solution to the asymptotics (|57p . 

Figure 1 shows the first three eigenvalues as functions of a obtained in three different ways. The solid 
lines are numerically computed from shooting method, the dashed lines are obtained using the expressions 
in (|2Tj) . and the dash-dot lines (almost visually indistinguishable from solid lines) are obtained by (l20l) . It is 
seen that accuracy of numerical solutions compare with approximate analytical results is very high. 
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7 Conclusions 

We summarize main results of this paper: 

• We found low-lying spectrum for a class of operators which appear in the linearization of the simplified 
critical Keller - Segel system around the one-parameter family of stationary solutions. These operators 
have one negative and one near zero eigenvalue and as a result - as discussed in the introduction - 
the blowup asymptotics will be governed by a two-parameter deformation of the static solution. The 
construction of such deformations and ensuing results are outlined in the introduction. 

• We constructed a rigorous and robust version of the method of matched asymptotics. We believe it can 
be used a large variety of linear differential operators arising in the linearization of nonlinear equations 
and hopefully can be extended to nonlinear ones as well (in this case perturbation series, (|26f and (|39[) . 
should be replaced by fixed point iterations). 

There are two main limitations of our set-up: we deal with radially symmetric solutions and with 
adiabatic approximation ignoring evolution of chemical concentration. Hence we conclude by emphasizing 
the desirability of two further extensions of our analysis by considering 

• non-radially symmetric initial conditions; 

• the full Keller - Segel model (without the adiabatic approximation). 
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A Appendix: Solutions of jCq(/) = 



In this appendix we derive, using well-known properties of the confluent hypergeometric functions, some 
properties of solutions of equation C (f> = which were used in Section [3] 

Proposition 7. If X/2a ^positive integer, then there are two independent solutions, <j)Q and cf>i, of the 
equation Cq4> = 0, satisfying the bounds (|36l) and (|37[) . and having the Wronskian (|38p and the expansion 

61. 



Proof. The equation C (j) = is the eigenvalue equation for the spherically symmetric harmonic oscillator in 
D = 4: 



= 0. 



(51) 
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Changing in this equation both dependent and independent variables, 

we obtain the Kummer's (or a confluent hypergeometric) differential equation [60] 

z X " + (2 - z) X ' + (J^ - lj X = 0. 
Assuming that X/2a ^positive integer, the latter equation has two linearly independent solutions 

Xo = U(l-^-,2,z), xi =1^(1-— ,2,*), 



(52) 



(53) 



(54) 



where U(a,b,z) and iFi(a,b, z) are the confluent hypergeometric functions of the second kind and the first 
kind, respectively [BD]- They are given by the following expressions (see the equations (13.1.2), (13.1.6) in 
[SO] and the equation (13) of the section 6.7.1 of [6Tj): 



1 F 1 (a,b,z) = ^2 



U(a,n + 1, z) — 



(a) r z r 

(_!)"+! 



n\T(a — n) 
(a) r z r 



(n + l) r r! 



i-Fi(a, n + 1, z)lnz 
{*(a + r) - *(1 + r) - *(1 + n + r)} 



(55) 



+ 



(n — 1)! ^— ^ (a — n) r z 



. a 



(1 — n) r r\ 



where (a)j — a(a + l)(a + 2) . . . (a + j — 1), (a)o = 1, n = 0, 1, 2, . . ., and the principal branch of In z is 
assumed to be chosen by setting — tt < argz < tt. 



The equations (|52[) and (|54[) give two linearly independent solutions 4>o and 0i of C (f> = on [0,oo) 

la la I 

My) = 1 f 1 (-^ + i,2,^)e^ 2 / 4; 



(56) 



where, using that <po and 4>\ are defined up to arbitrary constants, we added, for convenience, the factor 
r(— ^-). (Without that factor the equation (|44|) would have a factor l/r(— in all terms.) 

Asymptotic expansions of both confluent hypergeometric functions in (|54l) for ay 2 — > 00 are given by 
(see the equations (13.1.4) and (13.1.8) of Ref. [gOJ) 



2\ £-1 



1 + 



ay 



My) = 



r(i-£) V 2 



ay 



1 + ( — 

ay 



(57) 
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It follows from equations (1571) that out of these two functions only tfio decays at infinity (with correct 
asymptotic 4>q oc y^' 2 e~^ v ). 

The bounds (136[) and (|37[) are proven similarly, so we prove only (|36l) . Recall the definition of \o m 
(|54|) . ([55]) implies that for a constant C\, which might depend only on A/2a, there is a point zq > 1 s.t. 
|Xo(z)| < Ci-z^ -1 for all z > zq. Since xo depends only on z and X/2a, there is a constant C2, which might 
depend only on A/2a, s.t. |xo( z )| < C^z^ -1 for all z > 1. 

Now, ([551) implies that |xo( z )l < C^/z for all < z < 1 for some absolute constant C3. Combining the 
above inequalities gives the bound (|36l) . 

Now we compute the Wronskian, W(4>o,(f>i) :— <fiod y <f)i — dy<fio<f>i, of the two solutions found above. As 
usual, using the equation Cq(/) = 0, we derive the first order equation, W = — |W, for W(4>o, 0i). Solving 
this equation, we obtain W {<po , 4>\) — with a real constant C. To find this constant we compute W as 
y —> 00, using the equations ([57)1 and the fact T(z + 1) = zT(z) (on can also find C using the expansion 
d55j). This gives C = -f and (|38]>. 

Finally, we prove (fHT) for 0o- To this end we study the behaviour of the solution for y <§C 1/ai In 5 i, 
a <C 1, or equivalently, z <C 1/ln-. In the small z region, 

2 2a „ , . 
XA (^)=ln-- — +0(z) 

and e _ ^ z = 1 — ^z + O (z 2 ) . Computing the product of the small z expansions of %a and e~i , and replacing 
z with in the result gives the expression (|44|) . 

Now let y £ (0, Using the expansion (|55|) for n = 1 and ((52^) we obtain (|^)) . □ 
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